Tree topology diagnostics

Long-billed hermit song cultural evolution

Author
Published

October 10, 2024

Source code and data found at https://github.com/maRce10/lbh_cultural_evolution_revBayes

Load packages

Custom functions and global options

Code
knitr::opts_chunk$set(out.width = "100%", out.height = "100%", echo = TRUE) 

source('~/Dropbox/Projects/lbh_cultural_evolution_revBayes/source/custom_treespace.R')

theme_set(theme_classic(base_size = 10))

fill_color <- viridis::mako(10, alpha = 0.5)[7]

iterations <- 5000
chains <- 4
priors <- c(set_prior("normal(0, 1.5)", class = "Intercept"), set_prior("normal(0, 1.5)", class = "b"))

# to create several posterior predictive check plots out of a brms fit
custom_ppc1 <- function(fit, group = NULL, ndraws = 30) {
  
    ppc_dens <-
        pp_check(fit,
                 ndraws = ndraws,
                 type = 'dens_overlay_grouped',
                 group = group)  # shows dens_overlay plot by default
    pp_mean <-
        pp_check(
            fit,
            type = "stat_grouped",
            stat = "mean",
            group = group,
            ndraws = ndraws
        )
    pp_scat <-
        pp_check(fit,
                 type = "scatter_avg_grouped",
                 group = group,
                 ndraws = ndraws)
    pp_stat2 <- pp_check(fit, type = "stat_2d", ndraws = ndraws)
    pp_loo <- pp_check(fit, type = "loo_pit_qq", ndraws = ndraws)
    pp_error <-
        pp_check(fit,
                 type = "error_scatter_avg_grouped",
                 ndraws = ndraws,
                 group = group)
    plot_l <-
        list(ppc_dens, pp_mean, pp_scat, pp_error,  pp_stat2, pp_loo)
    
    plot_l[c(1, 3:length(plot_l))] <-
        lapply(plot_l[c(1, 3:length(plot_l))], function(x)
            x  + scale_color_viridis_d(
                begin = 0.1,
                end = 0.8,
                alpha = 0.5
            ) + scale_fill_viridis_d(
                begin = 0.1,
                end = 0.8,
                alpha = 0.5
            ) + theme_classic())     

    print(plot_grid(plotlist = plot_l, ncol = 2))
}

custom_ppc <- function(...) suppressMessages(custom_ppc1(...))

Read data

Code
rb.trees <- readRDS("./output/revbayes_output_in_single_R_object.RDS")

rb.trees <- rb.trees[grep("run_", names(rb.trees), invert = TRUE, value = TRUE)]

trees_diags <- data.frame(model = names(rb.trees))

# get lek name
trees_diags$lek <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 1)

# substitution model
trees_diags$subs <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 2)

# period
trees_diags$period <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 3)

# and which fossils were used
trees_diags$fossils <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 4)

trees_diags$align <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 5)

trees_diags$n_trees <- sapply(rb.trees, function(x) length(x$trees))

# Number of models by iterations and lek
kbl <-
    knitr::kable(as.matrix(table(trees_diags$lek, trees_diags$n_trees)), caption = "Number of models with a specific number of trees by lek")

kableExtra::kable_styling(kbl)

# Number of models by iterations and lek
kbl <-
    knitr::kable(as.matrix(table(trees_diags$lek, trees_diags$subs)), caption = "Number of models with a specific number of trees by lek")

kableExtra::kable_styling(kbl)

1 Topological space estimation

  • Pairwise topological distance between all trees for each lek
  • 100 trees for each model evenly spaced along the chain
  • Comparing based on shared song types (pairwise shared tips)
  • Some trees were removed if NAs were produced when comparing topologies (i.e. non-comparable topologies)

1.1 All fossils

Code
#selected leks
sel_leks <- c("SUR", "CCE", "HC1", "BR1", "TR1")

tree_names <- grep("early", names(rb.trees), invert = TRUE, value = TRUE)

pnts_lks <- lapply(sel_leks, function(i) {
  
  print(i)
  lek_trees <- grep(i, tree_names, value = TRUE)
  
  multi_list <- lapply(lek_trees, function(x) read.tree(file.path("./output/most_recent_revbayes_models/", x), keep.multi = TRUE))
  
  names(multi_list) <- sapply(basename(lek_trees), function(x) paste(strsplit(x, "_", fixed = TRUE)[[1]][1:5], collapse = "_"))
  
  pnts <- try(custom_treespace(multi_list, n.points = 100, method = "RF", cl = 10), silent = TRUE)

  if (!is(pnts, "try-error")) 
    return(pnts) else
      return(NULL)
})

names(pnts_lks) <- sel_leks

sapply(pnts_lks, is.null)

saveRDS(pnts_lks, file = "./output/topology/tree_distance_all_fossils_pairwise_shared.RDS")

1.2 Early fossils

Code
#selected leks
sel_leks <- c("SUR", "CCE", "HC1", "BR1", "TR1")
# sel_leks <- c("SUR", "TR1")

tree_names <- grep("early", names(rb.trees), value = TRUE)

pnts_lks <- lapply(sel_leks, function(i) {
  
  print(i)
  lek_trees <- grep(i, tree_names, value = TRUE)
  
  multi_list <- lapply(lek_trees, function(x) read.tree(file.path("./output/most_recent_revbayes_models/", x), keep.multi = TRUE))
  
  names(multi_list) <- sapply(basename(lek_trees), function(x) paste(strsplit(x, "_", fixed = TRUE)[[1]][1:5], collapse = "_"))
  
  pnts <- try(custom_treespace(multi_list, n.points = 100, method = "RF", cl = 10), silent = TRUE)

  if (!is(pnts, "try-error")) 
    return(pnts) else
      return(NULL)
})

names(pnts_lks) <- sel_leks

sapply(pnts_lks, is.null)

saveRDS(pnts_lks, file = "./output/topology/tree_distance_early_fossils_pairwise_shared.RDS")

2 Plot topological spaces

2.1 All fossils

2.1.1 New and old data sets

Code
pnts_lks <- readRDS("./output/topology/tree_distance_all_fossils_pairwise_shared.RDS")

ggs <- lapply(pnts_lks, function(X) {

    ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, ncol = 2) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") + 
      ggtitle(label = sprintf("Tree space for %d trees", nrow(X) / length(unique(X$chain)))) 
    })

ggs
$SUR


$CCE


$HC1


$BR1


$TR1

2.1.2 Old data set

Code
ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("old", X$chain), ]   
  
  if (nrow(X) > 0)
  ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))
    })

ggs
$SUR


$CCE


$HC1


$BR1
NULL

$TR1
NULL

2.1.3 New data set

Code
ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("new", X$chain), ]   
  
  ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %d trees", nrow(X) / length(unique(X$chain))))

    })

ggs
$SUR


$CCE


$HC1


$BR1


$TR1

2.1.3.1 Clumping molecular clocks

2.1.3.1.1 Old data set
Code
ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("old", X$chain), ]   
  X$clock <- ifelse(grepl("global", X$chain), "Global", "Uexp")
  X$chain <- gsub("_global|_Uexp", "", X$chain)
  
    if (nrow(X) > 0)
    ggplot(data = X, aes(x = x, y = y)) +
    geom_path(alpha = 0.1, size = 0.75, show.legend = FALSE) +
    geom_point(aes(shape = clock,  col = clock), size = 4) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_color_viridis_d(option = "G", begin = 0.2, end = 0.8, alpha = 0.75) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %d trees", nrow(X) / length(unique(X$chain))))

    })

ggs
$SUR


$CCE


$HC1


$BR1
NULL

$TR1
NULL
2.1.3.1.2 New data set
Code
ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("new", X$chain), ]   
  X$clock <- ifelse(grepl("global", X$chain), "Global", "Uexp")
  X$chain <- gsub("_global|_Uexp", "", X$chain)
  
  ggplot(data = X, aes(x = x, y = y)) +
    geom_path(alpha = 0.1, size = 0.75, show.legend = FALSE) +
    geom_point(aes(shape = clock,  col = clock), size = 4) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_color_viridis_d(option = "G", begin = 0.2, end = 0.8, alpha = 0.75) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %d trees", nrow(X) / length(unique(X$chain))))

    })

ggs
$SUR


$CCE


$HC1


$BR1


$TR1

2.2 Early fossils

2.2.1 New and old data sets

Code
pnts_lks <- readRDS("./output/topology/tree_distance_early_fossils_pairwise_shared.RDS")

ggs <- lapply(pnts_lks, function(X) {

    ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") + 
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain)))) 

    })

ggs
$SUR


$CCE


$HC1


$BR1


$TR1

2.2.2 New data set

Code
ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("new", X$chain), ]   
  
  ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))
  })

ggs
$SUR


$CCE


$HC1


$BR1


$TR1

2.2.3 Old data set

Code
ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("old", X$chain), ]   
  
  if (nrow(X) > 0)
  ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))
  })

ggs
$SUR


$CCE


$HC1


$BR1
NULL

$TR1
NULL

2.2.4 Clumping molecular clocks

2.2.4.1 Old data set

Code
ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("old", X$chain), ]   
  X$clock <- ifelse(grepl("global", X$chain), "Global", "Uexp")
  X$chain <- gsub("_global|_Uexp", "", X$chain)
  
    if (nrow(X) > 0)
    ggplot(data = X, aes(x = x, y = y)) +
    geom_path(alpha = 0.1, size = 0.75, show.legend = FALSE) +
    geom_point(aes(shape = clock,  col = clock), size = 4) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_color_viridis_d(option = "G", begin = 0.2, end = 0.8, alpha = 0.75) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))

    })

ggs
$SUR


$CCE


$HC1


$BR1
NULL

$TR1
NULL

2.2.4.2 New data set

Code
ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("new", X$chain), ]   
  X$clock <- ifelse(grepl("global", X$chain), "Global", "Uexp")
  X$chain <- gsub("_global|_Uexp", "", X$chain)
  
  ggplot(data = X, aes(x = x, y = y)) +
    geom_path(alpha = 0.1, size = 0.75, show.legend = FALSE) +
    geom_point(aes(shape = clock,  col = clock), size = 4) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_color_viridis_d(option = "G", begin = 0.2, end = 0.8, alpha = 0.75) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))

    })

ggs
$SUR


$CCE


$HC1


$BR1


$TR1

3 Topological space dissimilarity

3.1 Observations per lek and alignment

All fossils

Code
# similarity
out <- lapply(pnts_lks, function(x){
  
  Y <- space_similarity(formula = chain ~ x + y, data = x, pb = FALSE, type = "density.overlap")
  
  Y$mean.overlap <- Y$mean.overlap #-  mean(Y$mean.overlap)
  # Y$overlap <- Y$overlap / max(Y$overlap)
  
  return(Y)  
})

topo_similarity_all <- do.call(rbind, out)

topo_similarity_all2 <- topo_similarity_all[, c("group.1",  "group.2", "overlap.1in2")]
names(topo_similarity_all2) <- c("group.1", "group.2", "mean.overlap")

topo_similarity_all3 <- topo_similarity_all[, c("group.2",  "group.1", "overlap.2in1")]
names(topo_similarity_all3) <- c("group.2", "group.1", "mean.overlap")

topo_similarity_all <- rbind(topo_similarity_all2, topo_similarity_all3)

topo_similarity_all$alignment <- sapply(topo_similarity_all$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][1])

topo_similarity_all$data.set <- sapply(topo_similarity_all$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][2])

topo_similarity_all$fossils <- sapply(topo_similarity_all$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][3])

topo_similarity_all$models <- sapply(topo_similarity_all$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][4])


topo_similarity_all$lek <- substr(rownames(topo_similarity_all), 0 , 3)

tab <- as.data.frame(table(topo_similarity_all$lek, topo_similarity_all$alignment))
names(tab) <- c("lek", "model", "count")

tab
lek model count
BR1 all.equal 18
CCE all.equal 76
HC1 all.equal 76
SUR all.equal 76
TR1 all.equal 18
BR1 optimal 10
CCE optimal 44
HC1 optimal 44
SUR optimal 44
TR1 optimal 10
BR1 prank 2
CCE prank 12
HC1 prank 12
SUR prank 12
TR1 prank 2

Early fossils

Code
out <- lapply(pnts_lks, function(x){
  
  Y <- space_similarity(formula = chain ~ x + y, data = x, pb = FALSE, type = "mean.density.overlap")
  
  Y$mean.overlap <- Y$mean.overlap -  mean(Y$mean.overlap)
  # Y$overlap <- Y$overlap / max(Y$overlap)
  
  return(Y)  
})

topo_similarity_early <- do.call(rbind, out)

topo_similarity_early2 <- topo_similarity_early[, c("group.1",  "group.2", "overlap.2in1")]
names(topo_similarity_early2) <- c("group.2", "group.1", "mean.overlap")

topo_similarity_early3 <- topo_similarity_early[, c("group.1",  "group.2", "overlap.1in2")]
names(topo_similarity_early3) <- c("group.2", "group.1", "mean.overlap")

topo_similarity_early <- rbind(topo_similarity_early2, topo_similarity_early3)

topo_similarity_early$alignment <- sapply(topo_similarity_early$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][1])

topo_similarity_early$data.set <- sapply(topo_similarity_early$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][2])

topo_similarity_early$fossils <- sapply(topo_similarity_early$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][3])

topo_similarity_early$models <- sapply(topo_similarity_early$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][4])

topo_similarity_early$lek <- substr(rownames(topo_similarity_early), 0 , 3)

tab <- as.data.frame(table(topo_similarity_early$lek, topo_similarity_early$alignment))
names(tab) <- c("lek", "model", "count")

tab
lek model count
BR1 all.equal 2
CCE all.equal 12
HC1 all.equal 12
SUR all.equal 12
TR1 all.equal 2
BR1 optimal 10
CCE optimal 44
HC1 optimal 44
SUR optimal 44
TR1 optimal 10
BR1 prank 18
CCE prank 76
HC1 prank 76
SUR prank 76
TR1 prank 18

3.2 Plot dissimilarity by treatment

Code
topo_similarity_all$fossils <- "all"
topo_similarity_early$fossils <- "early"

shared_columns <- intersect(names(topo_similarity_all), names(topo_similarity_early))

topo_similarity <- rbind(topo_similarity_all[, shared_columns], topo_similarity_early[, shared_columns])

topo_similarity$dissimilarity <- 1 - topo_similarity$mean.overlap

# raincloud plot:
ggplot(topo_similarity,
       aes(x = alignment, y = dissimilarity, fill = alignment)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    # fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(# fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA) +
    # add justified jitter from the {gghalves} package
    gghalves::geom_half_point(
      # color = fill_color,
      # draw jitter on the left
      side = "l",
      # control range of jitter
      range_scale = .4,
      # add some transparency
      alpha = .5,
      transformation = ggplot2::position_jitter(height = 0)
    ) +
      scale_fill_viridis_d(option = "G",
                           end = 0.8,
                           alpha = 0.6) +
           scale_x_discrete(labels = c(
        "natural" = "Natural",
        "artificial" = "Artificial",
        "simulated" = "Simulated"
      )) +
      theme(legend.position = "none") +
      facet_wrap(~ data.set + models + fossils, ncol = 2) +
  labs(x = "Alignment", y = "Topologic space dissimilarity") 

3.3 Stats

3.3.1 Dissimilarity of the topological space

Code
# topo_similarity$dissimilarity[topo_similarity$dissimilarity <= 0] <- 0.001 

topo_similarity$alignment <- factor(topo_similarity$alignment, levels = c("all.equal", "optimal", "prank"))

topo_similarity$data.set <- factor(topo_similarity$data.set, levels = c("new", "old"))

topo_similarity$models <- factor(topo_similarity$models, levels = c("Uexp", "global"))

topo_similarity$fossils <- factor(topo_similarity$fossils, levels = c("early", "all"))

mod <- brm(
  formula = dissimilarity ~ alignment + data.set + models + fossils + (1 | lek), 
  data = topo_similarity,
  iter = iterations,
  family = gaussian(),
  chains = chains,
  cores = chains,
   control = list(adapt_delta = 0.99, max_treedepth = 15),
  file = "./data/processed/regression_results_topological_similiarity",
  file_refit = "always")

mod <- add_criterion(mod, criterion = c("loo"))

null_mod <- brm(
  formula = dissimilarity ~ 1, 
  data = topo_similarity,
  iter = iterations,
  family = gaussian(),
  chains = chains,
  cores = chains,
   control = list(adapt_delta = 0.99, max_treedepth = 15),
  file = "./data/processed/null_regression_results_topological_similiarity",
  file_refit = "always"
)

null_mod <- add_criterion(null_mod, criterion = c("loo"))

Model selection:

Code
mod <- readRDS("./data/processed/regression_results_topological_similiarity.rds")

null_mod <- readRDS("./data/processed/null_regression_results_topological_similiarity.rds")

loo_diff <- loo::loo_compare(mod, null_mod)

as.data.frame(loo_diff[,1:2], row.names = c("model", "null model"))
elpd_diff se_diff
model 0.00000 0.000000
null model -33.31634 8.978687
Code
extended_summary(fit = mod,
    n.posterior = 1000,
    fill =  mako(10)[7],
    trace.palette = viridis,
    remove.intercepts = TRUE,
    highlight = TRUE,
    trace = TRUE,
    return = FALSE, 
    print.name = FALSE
)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 Intercept-student_t(3, 0.9, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) dissimilarity ~ alignment + data.set + models + fossils + (1 | lek) 5000 4 1 2500 2 (2e-04%) 0 6259.166 5572.907 2016324370
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_alignmentoptimal 0.003 -0.057 0.063 1.001 6969.127 7503.027
b_alignmentprank 0.125 0.054 0.197 1 6259.166 6972.068
b_data.setold 0.188 0.138 0.237 1 8340.606 6680.098
b_modelsglobal 0.011 -0.035 0.056 1 8505.640 5572.907
b_fossilsall 0.090 0.029 0.149 1 6927.156 6884.685

Contrasts for alignment strategies:

Code
# contrasts
contrasts(
    fit = mod,
    sort.levels = c("prank", "optimal", "all.equal"),
    predictor = "alignment",
    n.posterior = 2000,
    level.sep = " VS ",
    html.table = TRUE,
    plot = TRUE,
    highlight = TRUE,
    fill = mako(10)[7]
)
Contrasts Estimate Est.Error l-95% CI u-95% CI
1 prank VS optimal 0.122 0.031 0.061 0.183
2 prank VS all.equal 0.125 0.037 0.054 0.197
3 optimal VS all.equal 0.003 0.031 -0.057 0.063

Code
custom_ppc(fit = mod, group = "alignment")

Session information

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggridges_0.5.6       cowplot_1.1.3        brmsish_1.0.0       
 [4] brms_2.22.0          Rcpp_1.0.13          PhenotypeSpace_0.1.0
 [7] ape_5.8              ggplot2_3.5.1        kableExtra_1.4.0    
[10] knitr_1.48           viridis_0.6.5        viridisLite_0.4.2   
[13] pbapply_1.7-2       

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2.1       rstudioapi_0.16.0      jsonlite_1.8.9        
  [4] magrittr_2.0.3         TH.data_1.1-2          estimability_1.5.1    
  [7] spatstat.utils_3.1-0   farver_2.1.2           rmarkdown_2.28        
 [10] vctrs_0.6.5            spatstat.explore_3.3-2 terra_1.7-78          
 [13] htmltools_0.5.8.1      curl_5.2.2             distributional_0.5.0  
 [16] tidybayes_3.0.7        raster_3.6-26          StanHeaders_2.32.10   
 [19] KernSmooth_2.23-24     htmlwidgets_1.6.4      plyr_1.8.9            
 [22] sandwich_3.1-0         emmeans_1.10.3         zoo_1.8-12            
 [25] igraph_2.0.3           lifecycle_1.0.4        pkgconfig_2.0.3       
 [28] Matrix_1.7-0           R6_2.5.1               fastmap_1.2.0         
 [31] digest_0.6.37          colorspace_2.1-1       CircStats_0.2-6       
 [34] tensor_1.5             vegan_2.6-8            labeling_0.4.3        
 [37] fansi_1.0.6            spatstat.sparse_3.1-0  adehabitatHR_0.4.22   
 [40] polyclip_1.10-7        abind_1.4-8            mgcv_1.9-1            
 [43] compiler_4.4.1         proxy_0.4-27           withr_3.0.1           
 [46] backports_1.5.0        inline_0.3.19          DBI_1.2.3             
 [49] QuickJSR_1.4.0         pkgbuild_1.4.4         MASS_7.3-61           
 [52] classInt_0.4-10        loo_2.8.0              permute_0.9-7         
 [55] tools_4.4.1            units_0.8-5            goftest_1.2-3         
 [58] quadprog_1.5-8         glue_1.8.0             nlme_3.1-165          
 [61] grid_4.4.1             sf_1.0-17              checkmate_2.3.2       
 [64] reshape2_1.4.4         cluster_2.1.6          ade4_1.7-22           
 [67] generics_0.1.3         gtable_0.3.5           spatstat.data_3.1-2   
 [70] class_7.3-22           tidyr_1.3.1            sp_2.1-4              
 [73] xml2_1.3.6             utf8_1.2.4             spatstat.geom_3.3-2   
 [76] pillar_1.9.0           ggdist_3.3.2           stringr_1.5.1         
 [79] posterior_1.6.0        splines_4.4.1          dplyr_1.1.4           
 [82] lattice_0.22-6         gghalves_0.1.4         survival_3.7-0        
 [85] deldir_2.0-4           tidyselect_1.2.1       arrayhelpers_1.1-0    
 [88] gridExtra_2.3          V8_4.4.2               svglite_2.1.3         
 [91] stats4_4.4.1           xfun_0.48              adehabitatMA_0.3.17   
 [94] bridgesampling_1.1-2   matrixStats_1.4.1      adehabitatLT_0.3.28   
 [97] rstan_2.32.6           stringi_1.8.4          yaml_2.3.10           
[100] boot_1.3-30            evaluate_1.0.0         codetools_0.2-20      
[103] tibble_3.2.1           cli_3.6.3              RcppParallel_5.1.9    
[106] xtable_1.8-4           systemfonts_1.1.0      munsell_0.5.1         
[109] spatstat.random_3.3-1  coda_0.19-4.1          svUnit_1.0.6          
[112] spatstat.univar_3.0-1  parallel_4.4.1         rstantools_2.4.0      
[115] bayesplot_1.11.1       Brobdingnag_1.2-9      phangorn_2.11.1       
[118] mvtnorm_1.3-1          scales_1.3.0           e1071_1.7-16          
[121] purrr_1.0.2            rlang_1.1.4            fastmatch_1.1-4       
[124] multcomp_1.4-25